{
 "metadata": {
  "name": "",
  "signature": "sha256:c5aa103d8ff8b16675145abcd90d602d7e982001758c23ed2c0a702709b2642b"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Metric Learning with the Shogun Machine Learning Toolbox"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "*By Fernando J. Iglesias Garcia (GitHub ID: [iglesias](https://github.com/iglesias)) as project report for GSoC 2013 ([project details](http://www.google-melange.com/gsoc/project/google/gsoc2013/iglesias/62013)).*"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook illustrates <a href=\"http://en.wikipedia.org/wiki/Statistical_classification\">classification</a> and <a href=\"http://en.wikipedia.org/wiki/Feature_selection\">feature selection</a> using <a href=\"http://en.wikipedia.org/wiki/Similarity_learning#Metric_learning\">metric learning</a> in Shogun. To overcome the limitations of <a href=\"http://en.wikipedia.org/wiki/K-nearest_neighbors_algorithm\">knn</a> with Euclidean distance as the distance measure, <a href=\"http://en.wikipedia.org/wiki/Large_margin_nearest_neighbor\">Large Margin Nearest Neighbour</a>(LMNN) is discussed. This is consolidated by applying LMNN over the metagenomics data set."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Building up the intuition to understand LMNN"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First of all, let us introduce LMNN through a simple example. For this purpose, we will be using the following two-dimensional toy data set:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy\n",
      "\n",
      "x = numpy.array([[0,0],[-1,0.1],[0.3,-0.05],[0.7,0.3],[-0.2,-0.6],[-0.15,-0.63],[-0.25,0.55],[-0.28,0.67]])\n",
      "y = numpy.array([0,0,0,0,1,1,2,2])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That is, there are eight feature vectors where each of them belongs to one out of three different classes (identified by either 0, 1, or 2). Let us have a look at this data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as pyplot\n",
      "\n",
      "%matplotlib inline\n",
      "\n",
      "def plot_data(features,labels,axis,alpha=1.0):\n",
      "    # separate features according to their class\n",
      "    X0,X1,X2 = features[labels==0], features[labels==1], features[labels==2]\n",
      "    \n",
      "    # class 0 data\n",
      "    axis.plot(X0[:,0], X0[:,1], 'o', color='green', markersize=12, alpha=alpha)\n",
      "    # class 1 data\n",
      "    axis.plot(X1[:,0], X1[:,1], 'o', color='red', markersize=12, alpha=alpha)\n",
      "    # class 2 data\n",
      "    axis.plot(X2[:,0], X2[:,1], 'o', color='blue', markersize=12, alpha=alpha)\n",
      "    \n",
      "    # set axes limits\n",
      "    axis.set_xlim(-1.5,1.5)\n",
      "    axis.set_ylim(-1.5,1.5)\n",
      "    axis.set_aspect('equal')\n",
      "    \n",
      "    axis.set_xlabel('x')\n",
      "    axis.set_ylabel('y')\n",
      "\n",
      "figure,axis = pyplot.subplots(1,1)\n",
      "plot_data(x,y,axis)\n",
      "axis.set_title('Toy data set')\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the figure above, we can see that two of the classes are represented by two points that are, for each of these classes, very close to each other. The third class, however, has four points that are close to each other with respect to the y-axis, but spread along the x-axis.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If we were to apply kNN (*k-nearest neighbors*) in a data set like this, we would expect quite some errors using the standard Euclidean distance. This is due to the fact that the spread of the data is not similar amongst the feature dimensions. The following piece of code plots an ellipse on top of the data set. The ellipse in this case is in fact a circunference that helps to visualize how the Euclidean distance weights equally both feature dimensions.   "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def make_covariance_ellipse(covariance):\n",
      "    import matplotlib.patches as patches\n",
      "    import scipy.linalg       as linalg\n",
      "    \n",
      "    # the ellipse is centered at (0,0)\n",
      "    mean = numpy.array([0,0])\n",
      "    \n",
      "    # eigenvalue decomposition of the covariance matrix (w are eigenvalues and v eigenvectors),\n",
      "    # keeping only the real part\n",
      "    w,v = linalg.eigh(covariance)\n",
      "    # normalize the eigenvector corresponding to the largest eigenvalue\n",
      "    u = v[0]/linalg.norm(v[0])\n",
      "    # angle in degrees\n",
      "    angle = 180.0/numpy.pi*numpy.arctan(u[1]/u[0])\n",
      "    # fill Gaussian ellipse at 2 standard deviation\n",
      "    ellipse = patches.Ellipse(mean, 2*w[0]**0.5, 2*w[1]**0.5, 180+angle, color='orange', alpha=0.3)\n",
      "    \n",
      "    return ellipse\n",
      "\n",
      "# represent the Euclidean distance\n",
      "figure,axis = pyplot.subplots(1,1)\n",
      "plot_data(x,y,axis)\n",
      "ellipse = make_covariance_ellipse(numpy.eye(2))\n",
      "axis.add_artist(ellipse)\n",
      "axis.set_title('Euclidean distance')\n",
      "pyplot.show()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A possible workaround to improve the performance of kNN in a data set like this would be to input to the kNN routine a distance measure. For instance, in the example above a good distance measure would give more weight to the y-direction than to the x-direction to account for the large spread along the x-axis. Nonetheless, it would be nicer (and, in fact, much more useful in practice) if this distance could be learnt automatically from the data at hand. Actually, LMNN is based upon this principle: given a number of neighbours *k*, find the Mahalanobis distance measure which maximizes kNN accuracy (using the given value for *k*) in a training data set. As we usually do in machine learning, under the assumption that the training data is an accurate enough representation of the underlying process, the distance learnt will not only perform well in the training data, but also have good generalization properties.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, let us use the [LMNN class](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLMNN.html) implemented in Shogun to find the distance and plot its associated ellipse. If everything goes well, we will see that the new ellipse only overlaps with the data points of the green class."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First, we need to wrap the data into Shogun's feature and label objects:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import RealFeatures, MulticlassLabels\n",
      "\n",
      "features = RealFeatures(x.T)\n",
      "labels   = MulticlassLabels(y.astype(numpy.float64))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Secondly, perform LMNN training:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import LMNN\n",
      "\n",
      "# number of target neighbours per example\n",
      "k = 1\n",
      "\n",
      "lmnn = LMNN(features,labels,k)\n",
      "# set an initial transform as a start point of the optimization\n",
      "init_transform = numpy.eye(2)\n",
      "lmnn.set_maxiter(2000)\n",
      "lmnn.train(init_transform)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "LMNN is an iterative algorithm. The argument given to [`train`](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLMNN.html#ab1b8bbdb8390415ac3ae7dc655cb512d) represents the initial state of the solution. By default, if no argument is given, then LMNN uses [PCA](http://en.wikipedia.org/wiki/Principal_component_analysis) to obtain this initial value."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, we retrieve the distance measure learnt by LMNN during training and visualize it together with the data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# get the linear transform from LMNN\n",
      "L = lmnn.get_linear_transform()\n",
      "# square the linear transform to obtain the Mahalanobis distance matrix\n",
      "M = numpy.matrix(numpy.dot(L.T,L))\n",
      "\n",
      "# represent the distance given by LMNN\n",
      "figure,axis = pyplot.subplots(1,1)\n",
      "plot_data(x,y,axis)\n",
      "ellipse = make_covariance_ellipse(M.I)\n",
      "axis.add_artist(ellipse)\n",
      "axis.set_title('LMNN distance')\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Beyond the main idea"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "LMNN is one of the so-called linear metric learning methods. What this means is that we can understand LMNN's output in two different ways: on the one hand, as a distance measure, this was explained above; on the other hand, as a linear transformation of the input data. Like any other linear transformation, LMNN's output can be written as a matrix, that we will call $L$. In other words, if the input data is represented by the matrix $X$, then LMNN can be understood as the data transformation expressed by $X'=L X$. We use the convention that each column is a feature vector; thus, the number of rows of $X$ is equal to the input dimension of the data, and the number of columns is equal to the number of vectors.\n",
      "\n",
      "So far, so good. But, if the output of the same method can be interpreted in two different ways, then there must be a relation between them! And that is precisely the case! As mentioned above, the ellipses that were plotted in the previous section represent a distance measure. This distance measure can be thought of as a matrix $M$, being the distance between two vectors $\\vec{x_i}$ and $\\vec{x_j}$ equal to $d(\\vec{x_i},\\vec{x_j})=(\\vec{x_i}-\\vec{x_j})^T M (\\vec{x_i}-\\vec{x_j})$. In general, this type of matrices are known as *Mahalanobis* matrices. In LMNN, the matrix $M$ is precisely the 'square' of the linear transformation $L$, i.e. $M=L^T L$. Note that a direct consequence of this is that $M$ is guaranteed to be positive semi-definite (PSD), and therefore define a valid metric.\n",
      "\n",
      "This distance measure/linear transform duality in LMNN has its own advantages. An important one is that the optimization problem can go back and forth between the $L$ and the $M$ representations, giving raise to a very efficient solution."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us now visualize LMNN using the linear transform interpretation. In the following figure we have taken our original toy data, transform it using $L$ and plot both the before and after versions of the data together."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# project original data using L\n",
      "lx = numpy.dot(L,x.T)\n",
      "\n",
      "# represent the data in the projected space\n",
      "figure,axis = pyplot.subplots(1,1)\n",
      "plot_data(lx.T,y,axis)\n",
      "plot_data(x,y,axis,0.3)\n",
      "ellipse = make_covariance_ellipse(numpy.eye(2))\n",
      "axis.add_artist(ellipse)\n",
      "axis.set_title('LMNN\\'s linear transform')\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the figure above, the transparent points represent the original data and are shown to ease the visualization of the LMNN transformation. Note also that the ellipse plotted is the one corresponding to the common Euclidean distance. This is actually an important consideration: if we think of LMNN as a linear transformation, the distance considered in the projected space is the Euclidean distance, and no any Mahalanobis distance given by M. To sum up, we can think of LMNN as a linear transform of the input space, or as method to obtain a distance measure to be used in the input space. It is an error to apply **both** the projection **and** the learnt Mahalanobis distance.  "
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Neighbourhood graphs"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "An alternative way to visualize the effect of using the distance found by LMNN together with kNN consists of using neighbourhood graphs. Despite the fancy name, these are actually pretty simple. The idea is just to construct a graph in the Euclidean space, where the points in the data set are the nodes of the graph, and a directed edge from one point to another denotes that the destination node is the 1-nearest neighbour of the origin node. Of course, it is also possible to work with neighbourhood graphs where $k \\gt 1$. Here we have taken the simplification of $k = 1$ so that the forthcoming plots are not too cluttered."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us define a data set for which the Euclidean distance performs considerably bad. In this data set there are several levels or layers in the y-direction. Each layer is populated by points that belong to the same class spread along the x-direction. The layers are close to each other in pairs, whereas the spread along x is larger. Let us define a function to generate such a data set and have a look at it."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy\n",
      "import matplotlib.pyplot as pyplot\n",
      "%matplotlib inline\n",
      "\n",
      "def sandwich_data():\n",
      "    from numpy.random import normal\n",
      "    \n",
      "    # number of distinct classes\n",
      "    num_classes = 6\n",
      "    # number of points per class\n",
      "    num_points = 9\n",
      "    # distance between layers, the points of each class are in a layer\n",
      "    dist = 0.7\n",
      "    \n",
      "    # memory pre-allocation\n",
      "    x = numpy.zeros((num_classes*num_points, 2))\n",
      "    y = numpy.zeros(num_classes*num_points)\n",
      "    \n",
      "    for i,j in zip(xrange(num_classes), xrange(-num_classes//2, num_classes//2 + 1)):\n",
      "        for k,l in zip(xrange(num_points), xrange(-num_points//2, num_points//2 + 1)):\n",
      "            x[i*num_points + k, :] = numpy.array([normal(l, 0.1), normal(dist*j, 0.1)])\n",
      "            \n",
      "        y[i*num_points:i*num_points + num_points] = i\n",
      "        \n",
      "    return x,y\n",
      "    \n",
      "\n",
      "def plot_sandwich_data(x, y, axis=pyplot, cols=['r', 'b', 'g', 'm', 'k', 'y']):\n",
      "    for idx,val in enumerate(numpy.unique(y)):\n",
      "        xi = x[y==val]\n",
      "        axis.scatter(xi[:,0], xi[:,1], s=50, facecolors='none', edgecolors=cols[idx])\n",
      "\n",
      "x, y = sandwich_data()\n",
      "figure, axis = pyplot.subplots(1, 1, figsize=(5,5))\n",
      "plot_sandwich_data(x, y, axis)\n",
      "\n",
      "axis.set_aspect('equal')\n",
      "axis.set_title('\"Sandwich\" toy data set')\n",
      "axis.set_xlabel('x')\n",
      "axis.set_ylabel('y')\n",
      "\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let the fun begin now! In the following block of code, we create an instance of a kNN classifier, compute the nearest neighbours using the Euclidean distance and, afterwards, using the distance computed by LMNN. The data set in the space result of the linear transformation given by LMNN is also shown."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import KNN, EuclideanDistance, LMNN, RealFeatures, MulticlassLabels\n",
      "\n",
      "def plot_neighborhood_graph(x, nn, axis=pyplot, cols=['r', 'b', 'g', 'm', 'k', 'y']):\n",
      "\tfor i in xrange(x.shape[0]):\n",
      "\t\txs = [x[i,0], x[nn[1,i], 0]]\n",
      "\t\tys = [x[i,1], x[nn[1,i], 1]]\n",
      "\t\taxis.plot(xs, ys, cols[int(y[i])])\n",
      "\n",
      "features = RealFeatures(x.T)\n",
      "labels = MulticlassLabels(y)\n",
      "\n",
      "fig, axes = pyplot.subplots(1, 3, figsize=(15, 10))\n",
      "\n",
      "# use k = 2 instead of 1 because otherwise the method nearest_neighbors just returns the same\n",
      "# points as their own 1-nearest neighbours\n",
      "k = 2\n",
      "\n",
      "knn = KNN(k, EuclideanDistance(features, features), labels)\n",
      "\n",
      "plot_sandwich_data(x, y, axes[0])\n",
      "plot_neighborhood_graph(x, knn.nearest_neighbors(), axes[0])\n",
      "axes[0].set_title('Euclidean neighbourhood in the input space')\n",
      "\n",
      "lmnn = LMNN(features, labels, k)\n",
      "# set a large number of iterations. The data set is small so it does not cost a lot, and this way\n",
      "# we ensure a robust solution\n",
      "lmnn.set_maxiter(3000)\n",
      "lmnn.train()\n",
      "knn.set_distance(lmnn.get_distance())\n",
      "\n",
      "plot_sandwich_data(x, y, axes[1])\n",
      "plot_neighborhood_graph(x, knn.nearest_neighbors(), axes[1])\n",
      "axes[1].set_title('LMNN neighbourhood in the input space')\n",
      "\n",
      "# plot features in the transformed space, with the neighbourhood graph computed using the Euclidean distance\n",
      "L = lmnn.get_linear_transform()\n",
      "xl = numpy.dot(x, L.T)\n",
      "features = RealFeatures(xl.T)\n",
      "knn.set_distance(EuclideanDistance(features, features))\n",
      "\n",
      "plot_sandwich_data(xl, y, axes[2])\n",
      "plot_neighborhood_graph(xl, knn.nearest_neighbors(), axes[2])\n",
      "axes[2].set_ylim(-3, 2.5)\n",
      "axes[2].set_title('Euclidean neighbourhood in the transformed space')\n",
      "\n",
      "[axes[i].set_xlabel('x') for i in xrange(len(axes))]\n",
      "[axes[i].set_ylabel('y') for i in xrange(len(axes))]\n",
      "[axes[i].set_aspect('equal') for i in xrange(len(axes))]\n",
      "\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice how all the lines that go across the different layers in the left hand side figure have disappeared in the figure in the middle. Indeed, LMNN did a pretty good job here. The figure in the right hand side shows the disposition of the points in the transformed space; from which the neighbourhoods in the middle figure should be clear. In any case, this toy example is just an illustration to give an idea of the power of LMNN. In the next section we will see how after applying a couple methods for feature normalization (e.g. scaling, whitening) the Euclidean distance is not so sensitive against different feature scales."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Real data sets"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Feature selection in metagenomics"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Metagenomics is a modern field in charge of the study of the DNA of microorganisms. The data set we have chosen for this section contains information about three different types of apes; in particular, gorillas, chimpanzees, and bonobos. Taking an approach based on metagenomics, the main idea is to study the DNA of the microorganisms (e.g. bacteria) which live inside the body of the apes. Owing to the many chemical reactions produced by these microorganisms, it is not only the DNA of the host itself important when studying, for instance, sickness or health, but also the DNA of the microorganisms inhabitants."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First of all, let us load the ape data set. This data set contains features taken from the bacteria inhabitant in the gut of the apes."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import CSVFile, RealFeatures, MulticlassLabels\n",
      "\n",
      "ape_features = RealFeatures(CSVFile('../../../data/multiclass/fm_ape_gut.dat'))\n",
      "ape_labels = MulticlassLabels(CSVFile('../../../data/multiclass/label_ape_gut.dat'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It is of course important to have a good insight of the data we are dealing with. For instance, how many examples and different features do we have?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print('Number of examples = %d, number of features = %d.' % (ape_features.get_num_vectors(), ape_features.get_num_features()))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So, 1472 features! Those are quite many features indeed. In other words, the feature vectors at hand lie on a 1472-dimensional space. We cannot visualize in the input feature space how the feature vectors look like. However, in order to gain a little bit more of understanding of the data, we can apply dimension reduction, embed the feature vectors in a two-dimensional space, and plot the vectors in the embedded space. To this end, we are going to use one of the [many](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CEmbeddingConverter.html) methods for dimension reduction included in Shogun. In this case, we are using t-distributed stochastic neighbour embedding (or [t-dsne](http://jmlr.org/papers/v9/vandermaaten08a.html)). This method is particularly suited to produce low-dimensional embeddings (two or three dimensions) that are straightforward to visualize."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def visualize_tdsne(features, labels):\n",
      "    from modshogun import TDistributedStochasticNeighborEmbedding\n",
      "    \n",
      "    converter = TDistributedStochasticNeighborEmbedding()\n",
      "    converter.set_target_dim(2)\n",
      "    converter.set_perplexity(25)\n",
      "    \n",
      "    embedding = converter.embed(features)\n",
      "    \n",
      "    import matplotlib.pyplot as pyplot\n",
      "    % matplotlib inline\n",
      "    \n",
      "    x = embedding.get_feature_matrix()\n",
      "    y = labels.get_labels()\n",
      "    \n",
      "    pyplot.scatter(x[0, y==0], x[1, y==0], color='green')\n",
      "    pyplot.scatter(x[0, y==1], x[1, y==1], color='red')\n",
      "    pyplot.scatter(x[0, y==2], x[1, y==2], color='blue')\n",
      "    pyplot.show()\n",
      "    \n",
      "visualize_tdsne(ape_features, ape_labels)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the figure above, the green points represent chimpanzees, the red ones bonobos, and the blue points gorillas. Providing the results in the figure, we can rapidly draw the conclusion that the three classes of apes are somewhat easy to discriminate in the data set since the classes are more or less well separated in two dimensions. Note that t-dsne use randomness in the embedding process. Thus, the figure result of the experiment in the previous block of code will be different after different executions. Feel free to play around and observe the results after different runs! After this, it should be clear that the bonobos form most of the times a very compact cluster, whereas the chimpanzee and gorillas clusters are more spread. Also, there tends to be a chimpanzee (a green point) closer to the gorillas' cluster. This is probably a outlier in the data set."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Even before applying LMNN to the ape gut data set, let us apply kNN classification and study how it performs using the typical Euclidean distance. Furthermore, since this data set is rather small in terms of number of examples, the kNN error above may vary considerably (I have observed variation of almost 20% a few times) across different runs. To get a robust estimate of how kNN performs in the data set, we will perform cross-validation using [Shogun's framework](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CCrossValidation.html) for evaluation. This will give us a reliable result regarding how well kNN performs in this data set."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import KNN, EuclideanDistance\n",
      "from modshogun import StratifiedCrossValidationSplitting, CrossValidation\n",
      "from modshogun import CrossValidationResult, MulticlassAccuracy\n",
      "\n",
      "# set up the classifier\n",
      "knn = KNN()\n",
      "knn.set_k(3)\n",
      "knn.set_distance(EuclideanDistance())\n",
      "\n",
      "# set up 5-fold cross-validation\n",
      "splitting = StratifiedCrossValidationSplitting(ape_labels, 5)\n",
      "# evaluation method\n",
      "evaluator = MulticlassAccuracy()\n",
      "cross_validation = CrossValidation(knn, ape_features, ape_labels, splitting, evaluator)\n",
      "# locking is not supported for kNN, deactivate it to avoid an inoffensive warning\n",
      "cross_validation.set_autolock(False)\n",
      "# number of experiments, the more we do, the less variance in the result\n",
      "num_runs = 200\n",
      "cross_validation.set_num_runs(num_runs)\n",
      "\n",
      "# perform cross-validation and print the result!\n",
      "result = cross_validation.evaluate()\n",
      "result = CrossValidationResult.obtain_from_generic(result)\n",
      "print('kNN mean accuracy in a total of %d runs is %.4f.' % (num_runs, result.mean))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, we can say that KNN performs actually pretty well in this data set. The average test classification error is less than between 2%. This error rate is already low and we should not really expect a significant improvement applying LMNN. This ought not be a surprise. Recall that the points in this data set have more than one thousand features and, as we saw before in the dimension reduction experiment, only two dimensions in an embedded space were enough to discern arguably well the chimpanzees, gorillas and bonobos."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that we have used stratified splitting for cross-validation. Stratified splitting divides the folds used during cross-validation so that the proportion of the classes in the initial data set is approximately maintained for each of the folds. This is particular useful in *skewed* data sets, where the number of examples among classes varies significantly."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Nonetheless, LMNN may still turn out to be very useful in a data set like this one. Making a small modification of the vanilla LMNN algorithm, we can enforce that the linear transform found by LMNN is diagonal. This means that LMNN can be used to weight each of the features and, once the training is performed, read from these weights which features are relevant to apply kNN and which ones are not. This is indeed a form of *feature selection*. Using Shogun, it is extremely easy to switch to this so-called *diagonal* mode for LMNN: just call the method [`set_diagonal(use_diagonal)`](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLMNN.html#ad2f03dbad3ad08ab76aecbb656a486e6) with `use_diagonal` set to `True`."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following experiment takes about five minutes until it is completed (using Shogun Release, i.e. compiled with optimizations enabled). This is mostly due to the high dimension of the data (1492 features) and the fact that, during training, LMNN has to compute many outer products of feature vectors, which is a computation whose time complexity is proportional to the square of the number of features. For the illustration purposes of this notebook, in the following cell we are just going to use a small subset of all the features so that the training finishes faster."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import LMNN\n",
      "import numpy\n",
      "\n",
      "# to make training faster, use a portion of the features\n",
      "fm = ape_features.get_feature_matrix()\n",
      "ape_features_subset = RealFeatures(fm[:150, :])\n",
      "\n",
      "# number of targer neighbours in LMNN, here we just use the same value that was used for KNN before\n",
      "k = 3\n",
      "lmnn = LMNN(ape_features_subset, ape_labels, k)\n",
      "lmnn.set_diagonal(True)\n",
      "lmnn.set_maxiter(1000)\n",
      "init_transform = numpy.eye(ape_features_subset.get_num_features())\n",
      "lmnn.train(init_transform)\n",
      "\n",
      "diagonal = numpy.diag(lmnn.get_linear_transform())\n",
      "print('%d out of %d elements are non-zero.' % (numpy.sum(diagonal != 0), diagonal.size))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So only 64 out of the 150 first features are important according to the result transform! The rest of them have been given a weight exactly equal to zero, even if all of the features were weighted equally with a value of one at the beginnning of the training. In fact, if all the 1472 features were used, only about 158 would have received a non-zero weight. Please, feel free to experiment using all the features!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It is a fair question to ask how did we know that the maximum number of iterations in this experiment should be around 1200 iterations. Well, the truth is that we know this only because we have run this experiment with this same data beforehand, and we know that after this number of iterations the algorithm has converged. This is not something nice, and the ideal case would be if one could completely forget about this parameter, so that LMNN uses as many iterations as it needs until it converges. Nevertheless, this is not practical at least because of two reasons:\n",
      "\n",
      "- If you are dealing with many examples or with very high dimensional feature vectors, you might not want to wait until the algorithm converges and have a look at what LMNN has found before it has completely converged.\n",
      "- As with any other algorithm based on gradient descent, the termination criteria can be tricky. Let us illustrate this further:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as pyplot\n",
      "%matplotlib inline\n",
      "\n",
      "statistics = lmnn.get_statistics()\n",
      "pyplot.plot(statistics.obj.get())\n",
      "pyplot.grid(True)\n",
      "pyplot.xlabel('Number of iterations')\n",
      "pyplot.ylabel('LMNN objective')\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Along approximately the first three hundred iterations, there is not much variation in the objective. In other words, the objective curve is pretty much flat. If we are not careful and use termination criteria that are not demanding enough, training could be stopped at this point. This would be wrong, and might have terrible results as the training had not clearly converged yet at that moment."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In order to avoid disastrous situations, in Shogun we have implemented LMNN with really demanding criteria for automatic termination of the training process. Albeit, it is possible to tune the termination criteria using the methods [`set_stepsize_threshold`](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLMNN.html#a76b6914cf9d1a53b0c9ecd828c7edbcb) and [`set_obj_threshold`](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLMNN.html#af78c7dd9ed2307c0d53e7383cdc01a24). These methods can be used to modify the lower bound required in the step size and the increment in the objective (relative to its absolute value), respectively, to stop training. Also, it is possible to set a hard upper bound on the number of iterations using [`set_maxiter`](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLMNN.html#afcf319806eb710a0d9535fbeddf93795) as we have done above. In case the internal termination criteria did not fire before the maximum number of iterations was reached, you will receive a warning message, similar to the one shown above. This is not a synonym that the training went wrong; but it is strongly recommended at this event to have a look at the objective plot as we have done in the previous block of code."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Multiclass classification"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In addition to feature selection, LMNN can be of course used for multiclass classification. I like to think about LMNN in multiclass classification as a way to empower kNN. That is, the idea is basically to apply kNN using the distance found by LMNN $-$ in contrast with using one of the other most common distances, such as the Euclidean one. To this end we will use the wine data set from the [UCI Machine Learning repository](http://archive.ics.uci.edu/ml/datasets/Wine \"Wine data set\")."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import CSVFile, RealFeatures, MulticlassLabels\n",
      "\n",
      "wine_features = RealFeatures(CSVFile('../../../data/uci/wine/fm_wine.dat'))\n",
      "wine_labels = MulticlassLabels(CSVFile('../../../data/uci/wine/label_wine.dat'))\n",
      "\n",
      "assert(wine_features.get_num_vectors() == wine_labels.get_num_labels())\n",
      "print('%d feature vectors with %d features from %d different classes.' % (wine_features.get_num_vectors(), \\\n",
      "       wine_features.get_num_features(), wine_labels.get_num_classes()))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First, let us evaluate the performance of kNN in this data set using the same cross-validation setting used in the previous section:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import KNN, EuclideanDistance\n",
      "from modshogun import StratifiedCrossValidationSplitting, CrossValidation\n",
      "from modshogun import CrossValidationResult, MulticlassAccuracy\n",
      "import numpy\n",
      "\n",
      "# kNN classifier\n",
      "k = 5\n",
      "knn = KNN()\n",
      "knn.set_k(k)\n",
      "knn.set_distance(EuclideanDistance())\n",
      "\n",
      "splitting = StratifiedCrossValidationSplitting(wine_labels, 5)\n",
      "evaluator = MulticlassAccuracy()\n",
      "cross_validation = CrossValidation(knn, wine_features, wine_labels, splitting, evaluator)\n",
      "cross_validation.set_autolock(False)\n",
      "num_runs = 200\n",
      "cross_validation.set_num_runs(num_runs)\n",
      "\n",
      "result = CrossValidationResult.obtain_from_generic(cross_validation.evaluate())\n",
      "euclidean_means = numpy.zeros(3)\n",
      "euclidean_means[0] = result.mean\n",
      "\n",
      "print('kNN accuracy with the Euclidean distance %.4f.' % result.mean)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Seconly, we will use LMNN to find a distance measure and use it with kNN:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import LMNN\n",
      "\n",
      "# train LMNN\n",
      "lmnn = LMNN(wine_features, wine_labels, k)\n",
      "lmnn.set_maxiter(1500)\n",
      "lmnn.train()\n",
      "\n",
      "# evaluate kNN using the distance learnt by LMNN\n",
      "knn.set_distance(lmnn.get_distance())\n",
      "\n",
      "result = CrossValidationResult.obtain_from_generic(cross_validation.evaluate())\n",
      "lmnn_means = numpy.zeros(3)\n",
      "lmnn_means[0] = result.mean\n",
      "\n",
      "print('kNN accuracy with the distance obtained by LMNN %.4f.' % result.mean)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The warning is fine in this case, we have made sure that the objective variation was really small after 1500 iterations. In any case, do not hesitate to check it yourself studying the objective plot as it was shown in the previous section."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As the results point out, LMNN really helps here to achieve better classification performance. However, this comparison is not entirely fair since the Euclidean distance is very sensitive to the scaling that different feature dimensions may have, whereas LMNN can adjust to this during training. Let us have a closer look to this fact. Next, we are going to retrieve the feature matrix and see what are the maxima and minima for every dimension."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print('minima = ' + str(numpy.min(wine_features, axis=1)))\n",
      "print('maxima = ' + str(numpy.max(wine_features, axis=1)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Examine the second and the last dimensions, for instance. The second dimension has values ranging from 0.74 to 5.8, while the values of the last dimension range from 278 to 1680. This will cause that the Euclidean distance works specially wrong in this data set. You can realize of this considering that the total distance between two points will almost certainly just take into account the contributions of the dimensions with largest range."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In order to produce a more fair comparison, we will rescale the data so that all the feature dimensions are within the interval [0,1]. Luckily, there is a [preprocessor](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CPreprocessor.html) class in Shogun that makes this straightforward."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import RescaleFeatures\n",
      "\n",
      "# preprocess features so that all of them vary within [0,1]\n",
      "preprocessor = RescaleFeatures()\n",
      "preprocessor.init(wine_features)\n",
      "wine_features.add_preprocessor(preprocessor)\n",
      "wine_features.apply_preprocessor()\n",
      "\n",
      "# sanity check\n",
      "assert(numpy.min(wine_features) >= 0.0 and numpy.max(wine_features) <= 1.0)\n",
      "\n",
      "# perform kNN classification after the feature rescaling\n",
      "knn.set_distance(EuclideanDistance())\n",
      "result = CrossValidationResult.obtain_from_generic(cross_validation.evaluate())\n",
      "euclidean_means[1] = result.mean\n",
      "\n",
      "print('kNN accuracy with the Euclidean distance after feature rescaling %.4f.' % result.mean)\n",
      "\n",
      "# train kNN in the new features and classify with kNN\n",
      "lmnn.train()\n",
      "knn.set_distance(lmnn.get_distance())\n",
      "result = CrossValidationResult.obtain_from_generic(cross_validation.evaluate())\n",
      "lmnn_means[1] = result.mean\n",
      "\n",
      "print('kNN accuracy with the distance obtained by LMNN after feature rescaling %.4f.' % result.mean)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Another different preprocessing that can be applied to the data is called *whitening*. Whitening, which is explained in an [article in wikipedia](http://en.wikipedia.org/wiki/Whitening_transformation \"Whitening transform\"), transforms the covariance matrix of the data into the identity matrix."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.linalg as linalg\n",
      "\n",
      "# shorthand for the feature matrix -- this makes a copy of the feature matrix\n",
      "data = wine_features.get_feature_matrix()\n",
      "# remove mean\n",
      "data = data.T\n",
      "data-= numpy.mean(data, axis=0)\n",
      "# compute the square of the covariance matrix and its inverse\n",
      "M = linalg.sqrtm(numpy.cov(data.T))\n",
      "# keep only the real part, although the imaginary that pops up in the sqrtm operation should be equal to zero\n",
      "N = linalg.inv(M).real\n",
      "# apply whitening transform\n",
      "white_data = numpy.dot(N, data.T)\n",
      "wine_white_features = RealFeatures(white_data)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The covariance matrices before and after the transformation can be compared to see that the covariance really becomes the identity matrix."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as pyplot\n",
      "%matplotlib inline\n",
      "\n",
      "fig, axarr = pyplot.subplots(1,2)\n",
      "axarr[0].matshow(numpy.cov(wine_features))\n",
      "axarr[1].matshow(numpy.cov(wine_white_features))\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, we evaluate again the performance obtained with kNN using the Euclidean distance and the distance found by LMNN using the whitened features."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "wine_features = wine_white_features\n",
      "\n",
      "# perform kNN classification after whitening\n",
      "knn.set_distance(EuclideanDistance())\n",
      "result = CrossValidationResult.obtain_from_generic(cross_validation.evaluate())\n",
      "euclidean_means[2] = result.mean\n",
      "\n",
      "print('kNN accuracy with the Euclidean distance after whitening %.4f.' % result.mean)\n",
      "\n",
      "# train kNN in the new features and classify with kNN\n",
      "lmnn.train()\n",
      "knn.set_distance(lmnn.get_distance())\n",
      "result = CrossValidationResult.obtain_from_generic(cross_validation.evaluate())\n",
      "lmnn_means[2] = result.mean\n",
      "\n",
      "print('kNN accuracy with the distance obtained by LMNN after whitening %.4f.' % result.mean)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As it can be seen, it did not really help to whiten the features in this data set with respect to only applying feature rescaling; the accuracy was already rather large after rescaling. In any case, it is good to know that this transformation exists, as it can become useful with other data sets, or before applying other machine learning algorithms."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us summarize the results obtained in this section with a bar chart grouping the accuracy results by distance (Euclidean or the one found by LMNN), and feature preprocessing:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "assert(euclidean_means.shape[0] == lmnn_means.shape[0])\n",
      "N = euclidean_means.shape[0]\n",
      "# the x locations for the groups\n",
      "ind = 0.5*numpy.arange(N)\n",
      "# bar width\n",
      "width = 0.15\n",
      "figure, axes = pyplot.subplots()\n",
      "figure.set_size_inches(6, 5)\n",
      "\n",
      "euclidean_rects = axes.bar(ind, euclidean_means, width, color='y')\n",
      "lmnn_rects = axes.bar(ind+width, lmnn_means, width, color='r')\n",
      "\n",
      "# attach information to chart\n",
      "axes.set_ylabel('Accuracies')\n",
      "axes.set_ylim(top=1.4)\n",
      "axes.set_title('kNN accuracy by distance and feature preprocessing')\n",
      "axes.set_xticks(ind+width)\n",
      "axes.set_xticklabels(('Raw', 'Rescaling', 'Whitening'))\n",
      "axes.legend(( euclidean_rects[0], lmnn_rects[0]), ('Euclidean', 'LMNN'), loc='upper right')\n",
      "\n",
      "def autolabel(rects):\n",
      "    # attach text labels to bars\n",
      "    for rect in rects:\n",
      "        height = rect.get_height()\n",
      "        axes.text(rect.get_x()+rect.get_width()/2., 1.05*height, '%.3f' % height,\n",
      "                  ha='center', va='bottom')\n",
      "        \n",
      "autolabel(euclidean_rects)\n",
      "autolabel(lmnn_rects)\n",
      "\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "References"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- Weinberger, K. Q., Saul, L. K. Distance Metric Learning for Large Margin Nearest Neighbor Classification. [(Link to paper in JMLR)](http://jmlr.org/papers/v10/weinberger09a.html)."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}
